Weibull minimum distribution (weibull_min)#

The Weibull distribution (SciPy: scipy.stats.weibull_min) is a flexible continuous distribution on \([0,\infty)\) that is widely used to model lifetimes / time-to-failure.

Its most important practical feature is that its hazard rate can be decreasing, constant, or increasing depending on the shape parameter — making it a go-to model in reliability engineering and survival analysis.


Learning goals#

  • Write down the PDF/CDF (and survival/hazard) and connect them to intuition.

  • Interpret the shape and scale parameters via the hazard rate.

  • Derive moments (mean/variance) using Gamma-function integrals.

  • Sample from the distribution using a NumPy-only inverse-CDF method.

  • Use scipy.stats.weibull_min for evaluation, sampling, and fitting.

import numpy as np

import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio

import scipy
from scipy import optimize, special
from scipy.stats import weibull_min as weibull_min_dist

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=5, suppress=True)
rng = np.random.default_rng(42)

# Record versions for reproducibility (useful when numerical details matter).
VERSIONS = {"numpy": np.__version__, "scipy": scipy.__version__, "plotly": plotly.__version__}

1) Title & Classification#

  • Name: weibull_min (Weibull minimum distribution; SciPy: scipy.stats.weibull_min)

  • Type: Continuous

  • Support (standard form): \(x \in [0,\infty)\)

  • Parameter space (standard form): shape \(c>0\)

  • SciPy location/scale: loc \in \mathbb{R}, scale > 0 with $\(X = \text{loc} + \text{scale}\,Y, \qquad Y \sim \mathrm{WeibullMin}(c).\)$

Unless stated otherwise, this notebook uses the standard form (loc=0, scale=1).

Note on naming: this is the usual Weibull distribution used for positive lifetimes. SciPy also provides weibull_max, which is a reflected version used in some extreme-value contexts.

2) Intuition & Motivation#

What it models#

The Weibull distribution is a workhorse model for positive durations and lifetimes. A central concept in reliability is the hazard rate (instantaneous failure rate)

\[ h(x) = \lim_{\Delta\downarrow 0}\frac{\mathbb{P}(x\le X < x+\Delta \mid X\ge x)}{\Delta}. \]

For a Weibull with shape \(c\) and scale \(\lambda\) (and loc=0), the hazard is a simple power law:

\[ h(x) = \frac{c}{\lambda}\left(\frac{x}{\lambda}\right)^{c-1}. \]

This gives a clean interpretation:

  • \(c<1\): decreasing hazard (“infant mortality” / early failures)

  • \(c=1\): constant hazard (memoryless exponential case)

  • \(c>1\): increasing hazard (“wear-out” / aging)

Typical real-world use cases#

  • Reliability / life testing: time-to-failure of components, fatigue life.

  • Survival analysis: parametric survival model when hazards are monotone.

  • Wind speed and hydrology: positive-valued environmental measurements.

  • Material strength: weakest-link arguments often motivate Weibull-like models.

Relations to other distributions#

  • Exponential: if \(c=1\), then \(X\sim \mathrm{Exp}(\text{scale}=\lambda)\).

  • Rayleigh: if \(c=2\) and \(\lambda=\sqrt{2}\,\sigma\), then \(X\) is Rayleigh(\(\sigma\)).

  • Gumbel (minimum) via log transform: if \(X\sim\mathrm{Weibull}(c,\lambda)\) and \(Y=\log X\), then $\(F_Y(y)=\mathbb{P}(Y\le y)=1-\exp\{-\exp(c(y-\log\lambda))\},\)$ which is a Gumbel-min (left-skewed extreme value) distribution.

  • Generative story from an exponential: if \(T\sim\mathrm{Exp}(1)\) then $\(X = \lambda\,T^{1/c} \sim \mathrm{Weibull}(c,\lambda).\)$ This directly yields an efficient sampler.

3) Formal Definition#

We use the common reliability parameterization with shape \(c>0\) and scale \(\lambda>0\) (SciPy’s scale). In the standard form, \(\lambda=1\).

Let \(z = \frac{x-\text{loc}}{\lambda}\).

PDF#

For \(x\ge \text{loc}\),

\[ f(x; c,\lambda,\text{loc}) = \frac{c}{\lambda}\,z^{c-1}\,\exp\{-z^c\}, \qquad z=\frac{x-\text{loc}}{\lambda}. \]

CDF#

For \(x\ge \text{loc}\),

\[ F(x; c,\lambda,\text{loc}) = 1 - \exp\{-z^c\}. \]

Equivalently, the survival function is

\[ S(x) = 1-F(x)=\exp\{-z^c\}. \]

Hazard rate#

Whenever \(S(x)>0\) (i.e., for finite \(x\)),

\[ h(x) = \frac{f(x)}{S(x)} = \frac{c}{\lambda} z^{c-1}. \]

Quantile function (PPF)#

For \(0<q<1\),

\[ F^{-1}(q) = \text{loc} + \lambda\,\bigl(-\log(1-q)\bigr)^{1/c}. \]
def weibull_min_pdf(x: np.ndarray, c: float) -> np.ndarray:
    """PDF of the standard WeibullMin(c) distribution (loc=0, scale=1)."""
    x = np.asarray(x, dtype=float)
    if c <= 0:
        return np.full_like(x, np.nan, dtype=float)

    out = np.zeros_like(x, dtype=float)
    mask = x >= 0
    xm = x[mask]

    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        out[mask] = c * np.power(xm, c - 1.0) * np.exp(-np.power(xm, c))

    return out


def weibull_min_logpdf(x: np.ndarray, c: float) -> np.ndarray:
    """Log-PDF of the standard WeibullMin(c) distribution (stable for tiny densities)."""
    x = np.asarray(x, dtype=float)
    if c <= 0:
        return np.full_like(x, np.nan, dtype=float)

    out = np.full_like(x, -np.inf, dtype=float)

    # Strictly positive region.
    mask_pos = x > 0
    xp = x[mask_pos]
    out[mask_pos] = np.log(c) + (c - 1.0) * np.log(xp) - np.power(xp, c)

    # Boundary at 0.
    mask_zero = x == 0
    if np.any(mask_zero):
        if np.isclose(c, 1.0):
            out[mask_zero] = 0.0  # pdf(0)=1
        elif c < 1.0:
            out[mask_zero] = np.inf  # pdf(0)=+inf
        else:
            out[mask_zero] = -np.inf  # pdf(0)=0

    return out


def weibull_min_cdf(x: np.ndarray, c: float) -> np.ndarray:
    """CDF of the standard WeibullMin(c) distribution."""
    x = np.asarray(x, dtype=float)
    if c <= 0:
        return np.full_like(x, np.nan, dtype=float)

    out = np.zeros_like(x, dtype=float)
    mask = x >= 0
    xm = x[mask]

    # Use expm1 for accuracy near x=0: 1 - exp(-t) = -expm1(-t).
    t = np.power(xm, c)
    out[mask] = -np.expm1(-t)

    return out


def weibull_min_sf(x: np.ndarray, c: float) -> np.ndarray:
    """Survival function S(x)=P(X>x) for the standard WeibullMin(c)."""
    x = np.asarray(x, dtype=float)
    if c <= 0:
        return np.full_like(x, np.nan, dtype=float)

    out = np.ones_like(x, dtype=float)
    mask = x >= 0
    xm = x[mask]
    out[mask] = np.exp(-np.power(xm, c))

    return out


def weibull_min_hazard(x: np.ndarray, c: float) -> np.ndarray:
    """Hazard rate h(x)=f(x)/S(x) for the standard WeibullMin(c)."""
    x = np.asarray(x, dtype=float)
    if c <= 0:
        return np.full_like(x, np.nan, dtype=float)

    out = np.zeros_like(x, dtype=float)
    mask = x >= 0
    xm = x[mask]

    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        out[mask] = c * np.power(xm, c - 1.0)

    return out


def weibull_min_ppf(q: np.ndarray, c: float) -> np.ndarray:
    """Quantile function (inverse CDF) of the standard WeibullMin(c)."""
    q = np.asarray(q, dtype=float)
    if c <= 0:
        return np.full_like(q, np.nan, dtype=float)

    out = np.full_like(q, np.nan, dtype=float)

    mask = (q >= 0) & (q <= 1)
    qm = q[mask]

    out[mask] = np.power(-np.log1p(-qm), 1.0 / c)
    out[q == 0] = 0.0
    out[q == 1] = np.inf

    return out


def weibull_min_pdf_loc_scale(x: np.ndarray, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """PDF with SciPy-style loc/scale using the standard-form implementation."""
    x = np.asarray(x, dtype=float)
    if scale <= 0:
        return np.full_like(x, np.nan, dtype=float)

    z = (x - loc) / scale
    out = np.zeros_like(x, dtype=float)
    mask = z >= 0
    out[mask] = weibull_min_pdf(z[mask], c) / scale
    return out


def weibull_min_cdf_loc_scale(x: np.ndarray, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """CDF with SciPy-style loc/scale using the standard-form implementation."""
    x = np.asarray(x, dtype=float)
    if scale <= 0:
        return np.full_like(x, np.nan, dtype=float)

    z = (x - loc) / scale
    out = np.zeros_like(x, dtype=float)
    mask = z >= 0
    out[mask] = weibull_min_cdf(z[mask], c)
    return out


def weibull_min_ppf_loc_scale(q: np.ndarray, c: float, *, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """PPF with SciPy-style loc/scale using the standard-form implementation."""
    if scale <= 0:
        q = np.asarray(q, dtype=float)
        return np.full_like(q, np.nan, dtype=float)

    return loc + scale * weibull_min_ppf(q, c)
# Sanity check: our formulas match SciPy (standard form and loc/scale).

c = 1.7
x = np.logspace(-4, 2, 25)
q = np.linspace(0.01, 0.99, 9)

rv_std = weibull_min_dist(c)

assert np.allclose(weibull_min_pdf(x, c), rv_std.pdf(x))
assert np.allclose(weibull_min_cdf(x, c), rv_std.cdf(x))
assert np.allclose(weibull_min_ppf(q, c), rv_std.ppf(q))

# loc/scale
loc, scale = -0.3, 2.5
rv_ls = weibull_min_dist(c, loc=loc, scale=scale)

assert np.allclose(weibull_min_pdf_loc_scale(x, c, loc=loc, scale=scale), rv_ls.pdf(x))
assert np.allclose(weibull_min_cdf_loc_scale(x, c, loc=loc, scale=scale), rv_ls.cdf(x))
assert np.allclose(weibull_min_ppf_loc_scale(q, c, loc=loc, scale=scale), rv_ls.ppf(q))

# hazard matches f/S in standard form (avoid survival underflow in the far tail)
x_haz = np.logspace(-4, 1, 30)
haz = weibull_min_pdf(x_haz, c) / weibull_min_sf(x_haz, c)
assert np.allclose(haz, weibull_min_hazard(x_haz, c))

4) Moments & Properties#

A convenient property of the Weibull distribution is that all positive moments exist and have a clean Gamma-function form.

Raw moments#

If \(X\sim\mathrm{Weibull}(c,\lambda)\) with loc=0, then for any \(r>-c\),

\[ \mathbb{E}[X^r] = \lambda^r\,\Gamma\!\left(1+\frac{r}{c}\right), \]

where \(\Gamma(\cdot)\) is the Gamma function.

Mean and variance#

Let \(g_k = \Gamma(1+k/c)\). Then

\[ \mathbb{E}[X] = \lambda\,g_1, \qquad \mathrm{Var}(X) = \lambda^2\,(g_2 - g_1^2). \]

Skewness and kurtosis#

Using raw moments and central-moment identities, the third and fourth central moments are

\[ \mu_3 = \lambda^3\,(g_3 - 3 g_1 g_2 + 2 g_1^3), \]
\[ \mu_4 = \lambda^4\,(g_4 - 4 g_1 g_3 + 6 g_1^2 g_2 - 3 g_1^4). \]

Skewness and excess kurtosis are

\[ \gamma_1 = \frac{\mu_3}{\sigma^3}, \qquad \gamma_2 = \frac{\mu_4}{\sigma^4} - 3. \]

MGF / characteristic function#

There is no simple elementary closed form for the moment generating function

\[M_X(t)=\mathbb{E}[e^{tX}],\]

but it can be written as a power series using the moments:

\[ M_X(t) = \sum_{n=0}^{\infty} \frac{t^n}{n!}\,\mathbb{E}[X^n] = \sum_{n=0}^{\infty}\frac{(t\lambda)^n}{n!}\,\Gamma\!\left(1+\frac{n}{c}\right), \]

with a radius of convergence that depends on \(c\):

  • \(c>1\): \(M_X(t)\) exists for all real \(t\) (the tail is lighter than exponential).

  • \(c=1\): \(M_X(t)\) exists for \(t < 1/\lambda\) (exponential case).

  • \(0<c<1\): \(M_X(t)\) diverges for every \(t>0\) (tail is heavier than exponential).

The characteristic function \(\varphi_X(t)=\mathbb{E}[e^{itX}]\) exists for all real \(t\) (bounded integrand).

Entropy#

The differential entropy (for loc=0) has a simple closed form:

\[ h(X) = 1 + \log\left(\frac{\lambda}{c}\right) + \gamma\,\left(1-\frac{1}{c}\right), \]

where \(\gamma\approx 0.57721\) is the Euler–Mascheroni constant.

EULER_GAMMA = float(-special.digamma(1.0))  # Euler–Mascheroni constant γ

def weibull_min_raw_moment(r: float, c: float, *, scale: float = 1.0) -> float:
    """E[X^r] for Weibull(c, scale) with loc=0."""
    if c <= 0 or scale <= 0:
        return float("nan")
    return float((scale**r) * special.gamma(1.0 + r / c))


def weibull_min_mean(c: float, *, scale: float = 1.0, loc: float = 0.0) -> float:
    if c <= 0 or scale <= 0:
        return float("nan")
    return float(loc + scale * special.gamma(1.0 + 1.0 / c))


def weibull_min_variance(c: float, *, scale: float = 1.0) -> float:
    if c <= 0 or scale <= 0:
        return float("nan")
    g1 = special.gamma(1.0 + 1.0 / c)
    g2 = special.gamma(1.0 + 2.0 / c)
    return float((scale**2) * (g2 - g1**2))


def weibull_min_skewness(c: float, *, scale: float = 1.0) -> float:
    if c <= 0 or scale <= 0:
        return float("nan")
    g1 = special.gamma(1.0 + 1.0 / c)
    g2 = special.gamma(1.0 + 2.0 / c)
    g3 = special.gamma(1.0 + 3.0 / c)

    mu3 = (scale**3) * (g3 - 3.0 * g1 * g2 + 2.0 * g1**3)
    var = (scale**2) * (g2 - g1**2)
    return float(mu3 / (var ** 1.5))


def weibull_min_excess_kurtosis(c: float, *, scale: float = 1.0) -> float:
    if c <= 0 or scale <= 0:
        return float("nan")
    g1 = special.gamma(1.0 + 1.0 / c)
    g2 = special.gamma(1.0 + 2.0 / c)
    g3 = special.gamma(1.0 + 3.0 / c)
    g4 = special.gamma(1.0 + 4.0 / c)

    mu4 = (scale**4) * (g4 - 4.0 * g1 * g3 + 6.0 * (g1**2) * g2 - 3.0 * g1**4)
    var = (scale**2) * (g2 - g1**2)
    return float(mu4 / (var**2) - 3.0)


def weibull_min_entropy(c: float, *, scale: float = 1.0) -> float:
    if c <= 0 or scale <= 0:
        return float("nan")
    return float(1.0 + np.log(scale / c) + EULER_GAMMA * (1.0 - 1.0 / c))


# Compare to SciPy.
c = 1.5
scale = 2.0
rv = weibull_min_dist(c, scale=scale)

mean_sp, var_sp, skew_sp, kurt_sp = rv.stats(moments="mvsk")
entropy_sp = rv.entropy()

mean_f = weibull_min_mean(c, scale=scale)
var_f = weibull_min_variance(c, scale=scale)
skew_f = weibull_min_skewness(c, scale=scale)
kurt_f = weibull_min_excess_kurtosis(c, scale=scale)
entropy_f = weibull_min_entropy(c, scale=scale)

print("SciPy vs formulas (c=1.5, scale=2.0):")
print("  mean   ", float(mean_sp), "|", mean_f)
print("  var    ", float(var_sp), "|", var_f)
print("  skew   ", float(skew_sp), "|", skew_f)
print("  kurt   ", float(kurt_sp), "|", kurt_f)
print("  entropy", float(entropy_sp), "|", entropy_f)

assert np.allclose([mean_sp, var_sp, skew_sp, kurt_sp, entropy_sp], [mean_f, var_f, skew_f, kurt_f, entropy_f])
SciPy vs formulas (c=1.5, scale=2.0):
  mean    1.805490585901867 | 1.805490585901867
  var     1.5027611392557279 | 1.5027611392557279
  skew    1.0719865728909583 | 1.0719865728909592
  kurt    1.3904035615957744 | 1.3904035615957788
  entropy 1.4800872940856251 | 1.4800872940856251

5) Parameter Interpretation#

SciPy’s weibull_min uses:

  • c as the shape parameter (often written \(k\) or \(\beta\))

  • scale as the scale parameter (often written \(\lambda\) or \(\eta\))

  • loc as a location shift

Shape c#

  • Controls the hazard rate behavior: $\(h(x)=\frac{c}{\lambda}\left(\frac{x}{\lambda}\right)^{c-1}.\)$

  • Also controls the shape near zero:

    • if \(c<1\), the density blows up at \(x=0\) (many very small values)

    • if \(c=1\), the density is finite at \(x=0\) (exponential)

    • if \(c>1\), the density is 0 at \(x=0\) and has a mode at \(x>0\)

Scale scale = \lambda#

  • Stretches the distribution horizontally: if \(Y\sim\mathrm{Weibull}(c,1)\) then \(X=\lambda Y\sim\mathrm{Weibull}(c,\lambda)\).

  • For lifetimes, \(\lambda\) is a characteristic life: \(F(\lambda)=1-e^{-1}\approx 0.632\).

Location loc#

  • Shifts support: support becomes \(x\ge \text{loc}\).

  • Useful for modeling a minimum lifetime (or a measurement offset).

# Shape effects: PDF and hazard for different c (standard form).

x_pdf = np.linspace(1e-4, 4.0, 800)
x_haz = np.logspace(-3, 1, 600)

c_values = [0.5, 1.0, 1.5, 3.0]

fig = go.Figure()
for c in c_values:
    y = weibull_min_pdf(x_pdf, c)
    fig.add_trace(go.Scatter(x=x_pdf, y=y, mode="lines", name=f"c={c}"))

fig.update_layout(
    title="Weibull_min PDF for different shapes (scale=1)",
    xaxis_title="x",
    yaxis_title="f(x; c)",
)
fig.show()

fig = go.Figure()
for c in c_values:
    y = weibull_min_hazard(x_haz, c)
    fig.add_trace(go.Scatter(x=x_haz, y=y, mode="lines", name=f"c={c}"))

fig.update_layout(
    title="Hazard rate h(x)=c x^{c-1} for different shapes (scale=1)",
    xaxis_title="x",
    yaxis_title="h(x)",
)
fig.update_xaxes(type="log")
fig.show()

# Scale effects at fixed c.
c = 2.0
scales = [0.5, 1.0, 2.0]

x = np.linspace(0, 6, 900)
fig = go.Figure()
for scale in scales:
    y = weibull_min_pdf_loc_scale(x, c, loc=0.0, scale=scale)
    fig.add_trace(go.Scatter(x=x, y=y, mode="lines", name=f"scale={scale}"))

fig.update_layout(
    title="Scale stretches the distribution (fixed c=2)",
    xaxis_title="x",
    yaxis_title="f(x; c, scale)",
)
fig.show()
# CDF shapes (standard form).

x = np.linspace(0, 6.0, 800)
c_values = [0.5, 1.0, 1.5, 3.0]

fig = go.Figure()
for c in c_values:
    fig.add_trace(go.Scatter(x=x, y=weibull_min_cdf(x, c), mode="lines", name=f"c={c}"))

fig.update_layout(
    title="Weibull_min CDF for different shapes (scale=1)",
    xaxis_title="x",
    yaxis_title="F(x; c)",
)
fig.show()

6) Derivations#

We derive moments and the likelihood in the loc=0 case. (Location just shifts \(X\) and does not change variance or shape.)

6.1 Expectation and general moments#

Start from the PDF with loc=0:

\[ f(x)=\frac{c}{\lambda}\left(\frac{x}{\lambda}\right)^{c-1}\exp\{-(x/\lambda)^c\},\qquad x\ge 0. \]

For \(r>-c\),

\[ \mathbb{E}[X^r] =\int_0^\infty x^r f(x)\,dx. \]

Use the substitution \(u=(x/\lambda)^c \Rightarrow x=\lambda u^{1/c}\) and \(dx = \lambda\,\frac{1}{c}\,u^{1/c - 1}\,du\).

After cancellation, the integral becomes

\[ \mathbb{E}[X^r] = \lambda^r \int_0^\infty u^{r/c} e^{-u}\,du = \lambda^r\,\Gamma\!\left(1+\frac{r}{c}\right). \]

Setting \(r=1\) and \(r=2\) yields mean and variance.

6.2 Variance#

Using \(\mathrm{Var}(X)=\mathbb{E}[X^2]-\mathbb{E}[X]^2\) and the moment formula:

\[ \mathrm{Var}(X) = \lambda^2\Bigl(\Gamma(1+2/c) - \Gamma(1+1/c)^2\Bigr). \]

6.3 Likelihood (i.i.d. sample)#

Let \(x_1,\dots,x_n\) be i.i.d. from a Weibull with parameters \((c,\lambda)\) and loc=0.

The likelihood is

\[ L(c,\lambda; x_{1:n}) = \prod_{i=1}^n \frac{c}{\lambda}\left(\frac{x_i}{\lambda}\right)^{c-1}\exp\{-(x_i/\lambda)^c\}. \]

The log-likelihood simplifies to

\[ \ell(c,\lambda) = n\log c + (c-1)\sum_{i=1}^n \log x_i - nc\log\lambda - \sum_{i=1}^n (x_i/\lambda)^c. \]

A useful fact: for fixed \(c\), the MLE of \(\lambda\) has a closed form:

\[ \hat\lambda(c)=\left(\frac{1}{n}\sum_{i=1}^n x_i^c\right)^{1/c}. \]

Plugging \(\hat\lambda(c)\) back into \(\ell\) yields a profile likelihood in \(c\); the resulting score equation has to be solved numerically.

def weibull_min_loglik(c: float, scale: float, x: np.ndarray) -> float:
    x = np.asarray(x, dtype=float)
    if c <= 0 or scale <= 0 or np.any(x <= 0):
        return -np.inf

    n = x.size
    return float(n * np.log(c) + (c - 1.0) * np.sum(np.log(x)) - n * c * np.log(scale) - np.sum((x / scale) ** c))


def weibull_min_scale_hat(c: float, x: np.ndarray) -> float:
    x = np.asarray(x, dtype=float)
    if c <= 0 or np.any(x <= 0):
        return float("nan")
    return float(np.mean(x**c) ** (1.0 / c))


def weibull_min_shape_score_profile(c: float, x: np.ndarray) -> float:
    """Score equation in c after profiling out scale (loc=0).

    Root of this function gives the MLE for c.
    """
    x = np.asarray(x, dtype=float)
    if c <= 0 or np.any(x <= 0):
        return float("nan")

    logx = np.log(x)
    x_c = x**c

    return float(1.0 / c + np.mean(logx) - np.sum(x_c * logx) / np.sum(x_c))


# MLE demo in the standard loc=0 form.
c_true = 1.7
scale_true = 2.5
x = weibull_min_dist(c_true, scale=scale_true).rvs(size=4_000, random_state=rng)

# Find a bracket where the profile score changes sign.
c_grid = np.linspace(0.15, 8.0, 400)
score_vals = np.array([weibull_min_shape_score_profile(c, x) for c in c_grid])

idx = np.where(np.sign(score_vals[:-1]) * np.sign(score_vals[1:]) < 0)[0]
if idx.size == 0:
    raise RuntimeError("Could not bracket the MLE root for c; try a wider grid.")

c_lo, c_hi = float(c_grid[idx[0]]), float(c_grid[idx[0] + 1])
sol = optimize.root_scalar(weibull_min_shape_score_profile, bracket=(c_lo, c_hi), args=(x,), method="brentq")

c_hat = float(sol.root)
scale_hat = weibull_min_scale_hat(c_hat, x)

print("True (c, scale):", (c_true, scale_true))
print("MLE  (c, scale):", (c_hat, scale_hat))

# Compare to SciPy's fit (fix loc=0).
c_hat_sp, loc_hat_sp, scale_hat_sp = weibull_min_dist.fit(x, floc=0)
print("SciPy fit (floc=0):", (float(c_hat_sp), float(scale_hat_sp)))

# Profile log-likelihood over c (scale profiled out).
c_grid = np.linspace(0.3, 5.0, 250)
ll_prof = np.array([weibull_min_loglik(c, weibull_min_scale_hat(c, x), x) for c in c_grid])

fig = go.Figure(go.Scatter(x=c_grid, y=ll_prof, mode="lines", name="profile loglik"))
fig.add_vline(x=c_true, line_dash="dash", line_color="green", annotation_text="true c")
fig.add_vline(x=c_hat, line_dash="dash", line_color="red", annotation_text="MLE c")
fig.update_layout(title="Profile log-likelihood for c (loc=0)", xaxis_title="c", yaxis_title="log-likelihood")
fig.show()
True (c, scale): (1.7, 2.5)
MLE  (c, scale): (1.681490527901736, 2.4991075083424734)
SciPy fit (floc=0): (1.6814619463263676, 2.499116178166809)

7) Sampling & Simulation#

NumPy-only algorithm (inverse transform)#

From the CDF in the loc=0 case:

\[ F(x)=1-\exp\{-(x/\lambda)^c\}. \]

Let \(U\sim\mathrm{Uniform}(0,1)\). Setting \(U=F(X)\) and solving for \(X\) gives

\[ X = \lambda\,\bigl(-\log(1-U)\bigr)^{1/c}. \]

This is an exact sampler (no rejection needed) and is typically the fastest way to generate Weibull samples.

def weibull_min_rvs_numpy(
    c: float,
    *,
    loc: float = 0.0,
    scale: float = 1.0,
    size: int | tuple[int, ...] = 1,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Sample from weibull_min using NumPy only (inverse-CDF sampler)."""
    if rng is None:
        rng = np.random.default_rng()

    if c <= 0 or scale <= 0:
        raise ValueError("Require c>0 and scale>0")

    u = rng.random(size=size)
    # -log1p(-u) is stable when u is very close to 1.
    x = scale * np.power(-np.log1p(-u), 1.0 / c)
    return loc + x

8) Visualization#

We’ll compare:

  • the theoretical PDF and CDF

  • Monte Carlo samples (NumPy-only sampler)

  • SciPy’s implementation

c = 1.3
scale = 2.0
n = 80_000

x_np = weibull_min_rvs_numpy(c, scale=scale, size=n, rng=rng)
x_sp = weibull_min_dist(c, scale=scale).rvs(size=n, random_state=rng)

# Histogram vs theoretical PDF
x_grid = np.linspace(0, np.quantile(x_np, 0.995), 500)
pdf_grid = weibull_min_pdf_loc_scale(x_grid, c, loc=0.0, scale=scale)

fig = px.histogram(
    x=x_np,
    nbins=140,
    histnorm="probability density",
    title="Monte Carlo histogram (NumPy-only) vs theoretical PDF",
    labels={"x": "x"},
)
fig.add_trace(go.Scatter(x=x_grid, y=pdf_grid, mode="lines", name="theoretical PDF"))
fig.show()

# Empirical CDF vs theoretical CDF
x_sorted = np.sort(x_np)
ecdf = np.arange(1, n + 1) / n
cdf_grid = weibull_min_cdf_loc_scale(x_grid, c, loc=0.0, scale=scale)

fig = go.Figure()
fig.add_trace(go.Scatter(x=x_sorted, y=ecdf, mode="lines", name="empirical CDF (NumPy-only)"))
fig.add_trace(go.Scatter(x=x_grid, y=cdf_grid, mode="lines", name="theoretical CDF"))
fig.update_layout(title="CDF: empirical vs theoretical", xaxis_title="x", yaxis_title="F(x)")
fig.show()

# Quick check: NumPy-only samples and SciPy samples should look similar.
from scipy.stats import ks_2samp

ks = ks_2samp(x_np, x_sp)
print("KS two-sample test (NumPy vs SciPy samples):")
print(ks)
KS two-sample test (NumPy vs SciPy samples):
KstestResult(statistic=0.004012499999999974, pvalue=0.5387311821802074, statistic_location=1.492943742605249, statistic_sign=-1)

9) SciPy Integration#

scipy.stats.weibull_min provides the standard distribution API:

  • weibull_min.pdf(x, c, loc=0, scale=1)

  • weibull_min.cdf(x, c, loc=0, scale=1)

  • weibull_min.rvs(c, loc=0, scale=1, size=..., random_state=...)

  • weibull_min.fit(data, ...) (MLE)

A common workflow is to freeze the distribution: rv = weibull_min(c, loc=..., scale=...), then call rv.pdf, rv.cdf, rv.rvs, etc.

c = 1.8
loc = 0.0
scale = 3.0

rv = weibull_min_dist(c, loc=loc, scale=scale)

x = np.array([0.2, 1.0, 3.0, 8.0])
print("pdf:", rv.pdf(x))
print("cdf:", rv.cdf(x))
print("sf :", rv.sf(x))

samples = rv.rvs(size=5, random_state=rng)
print("rvs:", samples)

# Fitting: estimate (c, scale) with loc fixed to 0.
true_c, true_scale = 1.4, 2.2
data = weibull_min_dist(true_c, scale=true_scale).rvs(size=5_000, random_state=rng)

c_hat, loc_hat, scale_hat = weibull_min_dist.fit(data, floc=0)
print("\nFit (fixed loc=0):")
print("  true (c, scale):", (true_c, true_scale))
print("  est  (c, scale):", (float(c_hat), float(scale_hat)))
pdf: [0.06823 0.21694 0.22073 0.00381]
cdf: [0.00761 0.12926 0.63212 0.9971 ]
sf : [0.99239 0.87074 0.36788 0.0029 ]
rvs: [2.36998 3.91892 3.90068 4.062   1.49376]

Fit (fixed loc=0):
  true (c, scale): (1.4, 2.2)
  est  (c, scale): (1.3901687026330698, 2.2085788623684626)

10) Statistical Use Cases#

Hypothesis testing (goodness-of-fit)#

If the parameters are specified in advance (not fit from the same sample), you can test whether data plausibly comes from a Weibull distribution using a goodness-of-fit test such as Kolmogorov–Smirnov (KS).

Caveat: if you estimate parameters from the data and then run KS on the same data, the usual KS p-values are no longer exact (use a parametric bootstrap or a corrected procedure).

Bayesian modeling#

There is no simple conjugate prior for \((c,\lambda)\) jointly, but there is a convenient conjugate update when shape \(c\) is known.

If \(X\sim\mathrm{Weibull}(c,\lambda)\) with loc=0, then \(Y=X^c\) has

\[ \mathbb{P}(Y\le y) = 1-\exp\{-y/\lambda^c\}, \]

so \(Y\sim\mathrm{Exp}(\text{rate}=\beta)\) with \(\beta = 1/\lambda^c\).

A Gamma prior on the rate \(\beta\) is conjugate.

Generative modeling#

Weibull distributions are commonly used as generative models for survival times and as components of mixture models (e.g., to model early-failure and wear-out subpopulations).

# Hypothesis testing example: KS test when parameters are known.

from scipy.stats import kstest

c = 1.2
scale = 2.0
x = weibull_min_dist(c, scale=scale).rvs(size=2_000, random_state=rng)

D, p_value = kstest(x, weibull_min_dist(c, scale=scale).cdf)
print("KS test against Weibull(c=1.2, scale=2.0):")
print("  D      =", D)
print("  p-value=", p_value)
KS test against Weibull(c=1.2, scale=2.0):
  D      = 0.03232217170754775
  p-value= 0.029962541380471608
# Bayesian modeling (conjugate update) when shape c is known.
# Model: X ~ Weibull(c, scale=lambda), loc=0.
# Transform: Y = X^c ~ Exp(rate=beta) with beta = 1 / lambda^c.
# Prior: beta ~ Gamma(alpha0, rate=r0).

from scipy.stats import gamma as gamma_dist

rng_local = np.random.default_rng(123)

c_known = 1.6
lambda_true = 2.5
x = weibull_min_dist(c_known, scale=lambda_true).rvs(size=800, random_state=rng_local)

y = x**c_known

# Conjugate Gamma prior on beta (rate parameterization).
alpha0 = 2.0
r0 = 1.0

alpha_post = alpha0 + y.size
r_post = r0 + np.sum(y)

# Draw posterior samples for beta, then transform to lambda.
beta_samps = gamma_dist(a=alpha_post, scale=1.0 / r_post).rvs(size=50_000, random_state=rng_local)
lambda_samps = (1.0 / beta_samps) ** (1.0 / c_known)

ci = np.quantile(lambda_samps, [0.05, 0.5, 0.95])
print("True lambda:", lambda_true)
print("Posterior lambda 90% CI + median:", ci)

fig = px.histogram(
    lambda_samps,
    nbins=120,
    histnorm="probability density",
    title="Posterior over scale (lambda) with known shape c",
    labels={"value": "lambda"},
)
fig.add_vline(x=lambda_true, line_dash="dash", line_color="green", annotation_text="true")
fig.show()
True lambda: 2.5
Posterior lambda 90% CI + median: [2.39838 2.48644 2.57815]
# Generative modeling example: early-failure vs wear-out mixture.

n = 60_000

# Early failures: c<1 (decreasing hazard), shorter characteristic life.
x_early = weibull_min_rvs_numpy(0.7, scale=0.8, size=n // 2, rng=rng)

# Wear-out: c>1 (increasing hazard), longer characteristic life.
x_wear = weibull_min_rvs_numpy(3.0, scale=2.0, size=n // 2, rng=rng)

x_mix = np.concatenate([x_early, x_wear])

fig = px.histogram(
    x_mix,
    nbins=160,
    histnorm="probability density",
    title="Mixture of Weibulls (early-failure + wear-out)",
    labels={"value": "time"},
)
fig.show()

print("Mixture summaries:")
print("  mean   =", float(x_mix.mean()))
print("  median =", float(np.median(x_mix)))
print("  90%    =", float(np.quantile(x_mix, 0.9)))
Mixture summaries:
  mean   = 1.4000915006124321
  median = 1.3065326608485597
  90%    = 2.6408964289216006

11) Pitfalls#

  • Invalid parameters: require c>0 and scale>0. With loc, support is \(x\ge \text{loc}\).

  • Boundary behavior at 0:

    • if \(c<1\), the PDF and hazard blow up at 0 (this is a feature of the model, not a bug).

    • plots that include exactly \(x=0\) can show infinities; start from a small \(\varepsilon>0\).

  • Numerical stability:

    • use logpdf when multiplying many densities or when probabilities are tiny;

    • for \(x\) near 0, compute CDF via -expm1(-t) (as done above);

    • for \(q\) near 1, use -log1p(-q) in the PPF (as done above).

  • Fitting with loc:

    • allowing loc to vary can lead to unstable fits or unintuitive parameter estimates;

    • in reliability, it’s common to fix loc=0 unless a physical minimum lifetime is justified.

  • Model misspecification:

    • Weibull enforces a monotone hazard; if the true hazard is bathtub-shaped (decrease then increase), consider mixtures or more flexible survival models.

12) Summary#

  • weibull_min is a continuous distribution on \([0,\infty)\) parameterized by shape \(c>0\) and (optionally) loc and scale.

  • Its CDF has the simple form \(F(x)=1-\exp\{-(x/\lambda)^c\}\), leading to an exact inverse-CDF sampler.

  • The hazard rate is \(h(x)=(c/\lambda)(x/\lambda)^{c-1}\): decreasing for \(c<1\), constant for \(c=1\), increasing for \(c>1\).

  • Moments are expressed with the Gamma function: \(\mathbb{E}[X^r]=\lambda^r\Gamma(1+r/c)\).

  • scipy.stats.weibull_min provides robust numerics for PDF/CDF/SF/PPF, sampling, and MLE fitting.

References#

  • Johnson, Kotz, and Balakrishnan. Continuous Univariate Distributions, Volume 1 (2nd ed.), Wiley, 1994.

  • Nelson. Applied Life Data Analysis, Wiley, 1982.

  • SciPy documentation: scipy.stats.weibull_min.